# setup
N = 40
# K = 2
from package.abn_mir_helper_functions import *
from package.abn_mir_plotting_functions import select_network
from package.bn_graph_methods import *
from package.plotting import get_colors, binary_states, bv_histogram
from random import SystemRandom as rand
from seq_pert_as_py import seq_pert_report, control_report
Problem specification:¶
Find a sequence of unit perturbations (p1, p2) such that for all cycles there exists an interval L such that a run-in from the perturbed state (i.e. p1 has been applied, L steps have occurred, and p2 has been applied) would terminate in cycle detection of the goal cycle.
Secondarily, let the Boolean network be a deterministic finite automaton that detects the goal cycle states and prevents further perturbation (in an applied model this could be something like cell death or apoptosis, for example).
Further,
Let there be a 5% chance of unit perturbation;
Alternatively,
Let there be a % chance of unit perturbation that is also a function of total change in Boolean states / |nodes| per step; proportional to total change / |nodes| or 1/total change / |nodes|
Solution method:¶
"It is better to be vaguely right than to be precisely wrong" - Carveth Read
Adjacency matrix from transitions matrix¶
After forming an adjacency matrix (with a 1 for each non-zero entry of the transitions matrix, and a 0 otherwise), by raising the matrix to successively higher powers, starting from 1, we could determine whether the goal states (cycle) are reachable from all other cycles, in some number of perturbations (the values of the non-zero entries after multiplication).
Similarly:
Algorithm:
For each perturbable node, form a matrix Apk, with elements ai,j, 1 if perturbation of node k causes a network in state i to transitions to state j, otherwise 0.
For all ordered pairs of matrices Apk:
For each row of Apk1 $*$ Apk2: Is the entry in the goal cycle column the largest value in the row? Count it as a good row.
Search for the maximum value of "good rows". (Later TODO: apply a different metric here)
If there is a non-zero entry in every value of the goal cycle column, then there does exist such a sequence of perturbations, separated by an interval L: 0 ≤ L ≤ longest cycle length - 1, such that applying p1, waiting L steps, and applying p2 will result in a goal state for every cycle. But, this interval may not be the same for each desired transition, and may not occur for every state of every cycle.
Set up a network for which to search for sequential perturbations effecting goal states¶
net = select_network(
num_nodes=N,
minimum_max_cycle_length=15,
maximum_max_cycle_length=35,
minimum_number_of_cycles=10,
maximum_number_of_cycles=30,
maximum_number_of_networks=1024,
iterations_limit=400,
)[
0
] # select net also returns all generated nets, used in bool_main for appendix figures
net.add_cycles(1000)
print(net)
net.compute_unit_perturbations_matrix(sort_selection=1, compute_over_t_u=False)
setup_colors = get_colors(
len(net.bn_collapsed_cycles) + 10, True
) # [[1, 1, 1, 1]] # [[0.3, 0.3, 0.3, 0.5]]
avg_color = [
sum(setup_colors[i][0] for i in range(len(setup_colors))) / len(setup_colors),
sum(setup_colors[i][1] for i in range(len(setup_colors))) / len(setup_colors),
sum(setup_colors[i][2] for i in range(len(setup_colors))) / len(setup_colors),
sum(setup_colors[i][3] for i in range(len(setup_colors))) / len(setup_colors) / 10,
]
setup_colors.append(avg_color)
BooleanNetwork: N=40. avg_k= 1.7. Run-ins: 1202. Observed cycles: 13. t_records: 1202. u_records: 0 Cycle lengths: [20, 20, 20, 20, 20, 20, 7, 4, 4, 4, 4, 4, 4]
Cycle states as binary numbers (left), goal cycle states as binary values (right; white: off, black: on)¶
Each state is taken to be a binary integer (ex: for N=2, a state 012 -> 1, a state 112 -> 3)
# goal states as a given cycle
goal_cycle_index = rand().randrange(0, len(net.bn_collapsed_cycles.cycle_records))
fig = binary_states(net, setup_colors, goal_cycle_index, plt)
plt.show()
Run-ins colored by cycle index without perturbation, without noise:¶
fig1, fig2 = control_report(
net=net,
cycle_colors=setup_colors,
total_steps=1000,
progress_div=1,
frequency_stacked=True,
with_noise=False,
# show_legend = True
)
plt.show()
Run-ins colored by cycle index without perturbation, with noise:¶
fig1, fig2 = control_report(
net=net,
cycle_colors=setup_colors,
total_steps=1000,
progress_div=1,
frequency_stacked=True,
with_noise=0.05,
show_legend=True,
)
plt.show()
Run-ins colored by cycle index without perturbation, with noise, BN = DFA:¶
fig1, fig2 = control_report(
net=net,
cycle_colors=setup_colors,
total_steps=1000,
progress_div=1,
frequency_stacked=True,
with_noise=0.05,
show_legend=True,
goal_cycle_index=goal_cycle_index,
goal_dfa = True
)
plt.show()
Applying algorithm to find a perturbation composed of two unit perturbations (=flipping the state of one bit)¶
Interval L then varies over finite state machine operation
matrices = [
np.zeros_like(net.cycles_unit_perturbations_transition_matrix) for _ in net.nodes
]
all_matrix = np.zeros_like(net.cycles_unit_perturbations_transition_matrix)
for record in net.cycles_unit_perturbations_records:
if record.end_index is not None:
matrices[record.perturbed_node_index][record.start_index][
record.end_index
] = 1 # 1 represents adjacency
all_matrix[record.start_index][
record.end_index
] = 1 # 1 represents adjacency # this one is for a secondary question
max_good_rows = 0
found_sequences_2 = []
# working_sequences_2 = [] # initial output -> thoughts: requiring that the perturbation not make the goal state change is important, but may be overly restrictive
best_sequence = None
for i in range(len(matrices)):
for j in range(len(matrices)):
matrix_product = matrices[i] @ matrices[j]
good_rows = 0
for row in matrix_product:
# if matrix_product[L][goal_cycle_index] == max(row):
# # if row[goal_cycle_index] > 0:
# good_rows += 1
good_rows += 1 # row[goal_cycle_index] / sum(row)
if good_rows == len(net.bn_collapsed_cycles) - 1:
found_sequences_2.append((i, j))
if good_rows > max_good_rows:
max_good_rows = good_rows
best_sequence = (i, j)
print(
f"{len(found_sequences_2)} sequences of perturbations unanimously effect goal state"
)
print(
f"results:\nsequence: {best_sequence}\nnumber of cycles for which sequence -> goal: {max_good_rows}"
)
0 sequences of perturbations unanimously effect goal state results: sequence: (0, 0) number of cycles for which sequence -> goal: 14
Applying algorithm to find a perturbation composed of three unit perturbations (=flipping the state of one bit)¶
Intervals L1, L2 then vary over finite state machine operation
matrices = [
np.zeros_like(net.cycles_unit_perturbations_transition_matrix) for _ in net.nodes
]
all_matrix = np.zeros_like(net.cycles_unit_perturbations_transition_matrix)
for record in net.cycles_unit_perturbations_records:
if record.end_index is not None:
matrices[record.perturbed_node_index][record.start_index][
record.end_index
] = 1 # 1 represents adjacency
all_matrix[record.start_index][
record.end_index
] = 1 # 1 represents adjacency # this one is for a secondary question
max_good_rows = 0
found_sequence = ""
found_sequences_3 = []
best_sequence_3 = None
# working_sequences_3 = []
for i in range(len(matrices)):
for j in range(len(matrices)):
m_i_j = matrices[i] @ matrices[j]
for k in range(len(matrices)):
matrix_product = m_i_j @ matrices[k]
good_rows = 0
# print(f"sum(row): {sum(row)}")
for row in matrix_product:
if row[goal_cycle_index] == max(row):
# if row[goal_cycle_index] > 0:
good_rows += 1 # row[goal_cycle_index] / sum(r for r in row)
# if good_rows == len(net.bn_collapsed_cycles) - 1:
# found_sequences_3.append((i, j, k))
if good_rows > max_good_rows:
max_good_rows = good_rows
found_sequence = (i, j, k)
print(
f"{len(found_sequences_3)} sequences of perturbations unanimously effect goal state"
)
print(
f"results:\nsequence: {found_sequence}\nnumber of cycles for which sequence -> goal: {max_good_rows}"
)
0 sequences of perturbations unanimously effect goal state results: sequence: (0, 2, 27) number of cycles for which sequence -> goal: 14
Progression of distribution of states from uniform initial conditions (1 for each cycle state), BN ≠DFA¶
With two unit perturbations, without noise, BN ≠DFA¶
fig1, fig2 = seq_pert_report(
p1=best_sequence[0],
p2=best_sequence[1],
goal_cycle_index=goal_cycle_index,
net=net,
cycle_colors=setup_colors,
total_steps=1000,
progress_div=1,
goal_bool=False,
frequency_stacked=True,
)
fig1.set_size_inches(10, 5)
fig2.set_size_inches(10, 5)
plt.show()
# color by terminal state label
goal cycle is in parentheses
With three unit perturbations, without noise, BN ≠DFA¶
fig1, fig2 = seq_pert_report(
p1=found_sequence[0],
p2=found_sequence[1],
goal_cycle_index=goal_cycle_index,
net=net,
cycle_colors=setup_colors,
total_steps=1000,
progress_div=1,
goal_bool=False,
p3=found_sequence[2],
frequency_stacked=True,
)
fig1.set_size_inches(10, 5)
fig2.set_size_inches(10, 5)
plt.show()
goal cycle is in parentheses
With noise:¶
With two unit perturbations, noise, BN ≠DFA:¶
fig1, fig2 = seq_pert_report(
p1=best_sequence[0],
p2=best_sequence[1],
goal_cycle_index=goal_cycle_index,
net=net,
cycle_colors=setup_colors,
total_steps=1000,
progress_div=1,
goal_bool=False,
with_noise=0.05,
frequency_stacked=True,
)
fig1.set_size_inches(10, 5)
fig2.set_size_inches(10, 5)
plt.show()
goal cycle is in parentheses
With three unit perturbations, noise, BN ≠DFA¶
fig1, fig2 = seq_pert_report(
p1=found_sequence[0],
p2=found_sequence[1],
goal_cycle_index=goal_cycle_index,
net=net,
cycle_colors=setup_colors,
total_steps=1000,
progress_div=1,
goal_bool=False,
p3=found_sequence[2],
with_noise=0.05,
frequency_stacked=True,
# show_legend = True,
)
fig1.set_size_inches(10, 5)
fig2.set_size_inches(10, 5)
plt.show()
goal cycle is in parentheses
Progression of distribution of states from uniform initial conditions (1 for each cycle state), BN = DFA¶
With two unit perturbations, without noise, BN = DFA¶
fig1, fig2 = seq_pert_report(
p1=best_sequence[0],
p2=best_sequence[1],
goal_cycle_index=goal_cycle_index,
net=net,
cycle_colors=setup_colors,
total_steps=1000,
progress_div=1,
goal_bool=True,
frequency_stacked=True,
)
fig1.set_size_inches(10, 5)
fig2.set_size_inches(10, 5)
plt.show()
goal cycle is in parentheses
With three unit perturbations, without noise, BN = DFA¶
fig1, fig2 = seq_pert_report(
p1=found_sequence[0],
p2=found_sequence[1],
goal_cycle_index=goal_cycle_index,
net=net,
cycle_colors=setup_colors,
total_steps=1000,
progress_div=1,
goal_bool=True,
p3=found_sequence[2],
frequency_stacked=True,
)
fig1.set_size_inches(10, 5)
fig2.set_size_inches(10, 5)
plt.show()
goal cycle is in parentheses
With noise:¶
With two unit perturbations, noise, BN = DFA¶
fig1, fig2 = seq_pert_report(
p1=best_sequence[0],
p2=best_sequence[1],
goal_cycle_index=goal_cycle_index,
net=net,
cycle_colors=setup_colors,
total_steps=1000,
progress_div=1,
goal_bool=True,
with_noise=0.05,
frequency_stacked=True,
)
fig1.set_size_inches(10, 5)
fig2.set_size_inches(10, 5)
plt.show()
goal cycle is in parentheses
With three unit perturbations, noise, BN = DFA¶
fig1, fig2 = seq_pert_report(
p1=found_sequence[0],
p2=found_sequence[1],
goal_cycle_index=goal_cycle_index,
net=net,
cycle_colors=setup_colors,
total_steps=1000,
progress_div=1,
goal_bool=True, # DFA
p3=found_sequence[2],
with_noise=0.05,
frequency_stacked=True,
)
fig1.set_size_inches(10, 5)
fig2.set_size_inches(10, 5)
plt.show()
goal cycle is in parentheses
With non-constant noise:¶
Step to step Hamming distance, or "Boolean velocity" (BV) noise (perturbations possible when more nodes are changing states):¶
If BV ≥ value: possible noise perturbation
Inverse BV noise (perturbation possible when more nodes are in constant states)¶
If BV < value: possible noise perturbation
# histogram of BV for this network
fig = bv_histogram(net, setup_colors, stacked=True, line_hist=True, norm_x=True)
plt.xlabel('Step to step Hamming distance ("Boolean velocity")')
plt.ylabel("number of observations")
plt.show()
BV noise:¶
2 unit perturbations; BN = DFA
fig1, fig2 = seq_pert_report(
p1=best_sequence[0],
p2=best_sequence[1],
goal_cycle_index=goal_cycle_index,
net=net,
cycle_colors=setup_colors,
total_steps=1000,
progress_div=1,
goal_bool=True, # DFA
with_noise="inverse bv",
frequency_stacked=True,
)
fig1.set_size_inches(10, 5)
fig2.set_size_inches(10, 5)
plt.show()
goal cycle is in parentheses
BV noise:¶
3 unit perturbations; BN = DFA
fig1, fig2 = seq_pert_report(
p1=found_sequence[0],
p2=found_sequence[1],
goal_cycle_index=goal_cycle_index,
net=net,
cycle_colors=setup_colors,
total_steps=1000,
progress_div=1,
goal_bool=True, # DFA
p3=found_sequence[2],
with_noise="bv",
frequency_stacked=True,
)
fig1.set_size_inches(10, 5)
fig2.set_size_inches(10, 5)
plt.show()
goal cycle is in parentheses
Inverse BV noise:¶
2 unit perturbations; BN = DFA
fig1, fig2 = seq_pert_report(
p1=best_sequence[0],
p2=best_sequence[1],
goal_cycle_index=goal_cycle_index,
net=net,
cycle_colors=setup_colors,
total_steps=1000,
progress_div=1,
goal_bool=True, # DFA
with_noise="inverse bv",
frequency_stacked=True,
)
fig1.set_size_inches(10, 5)
fig2.set_size_inches(10, 5)
plt.show()
goal cycle is in parentheses
Inverse BV noise:¶
3 unit perturbations; BN = DFA
fig1, fig2 = seq_pert_report(
p1=found_sequence[0],
p2=found_sequence[1],
goal_cycle_index=goal_cycle_index,
net=net,
cycle_colors=setup_colors,
total_steps=1000,
progress_div=1,
goal_bool=True, # DFA
p3=found_sequence[2],
with_noise="inverse bv",
frequency_stacked=True,
)
fig1.set_size_inches(10, 5)
fig2.set_size_inches(10, 5)
plt.show()
goal cycle is in parentheses